rm(list = ls(all = TRUE))
library(mvtnorm)
library(foreign)
library(blme)
d <- read.dta("leaders.dta")
load("rdmodels_main.Rdata")

## Functions
pfig <- function(mod, dat){
    out <- list(NULL)
    beta <- getME(mod, name = "beta") # fixed effects
    sigma <- as.matrix(vcov(mod)) # Variance covariance matrix
    beta.qb <- rmvnorm(1000, beta, sigma)
    dsim <- dat[dat$dep == 0, ]
    dsim$logdur <- log(dsim$duration + 1)
    uvals <- sort(unique(dsim$logdur))
    matsim <- cbind(rep(1, length(uvals)), uvals)
    colnames(matsim) <- c("intercept", "logduration")
    phat <- c()
    for(i in 1:nrow(beta.qb)){
        phat <- rbind(phat, t(exp(matsim %*% beta.qb[i, ]) / (1 + exp(matsim %*% beta.qb[i, ]))))
    }
    out[[1]] <- uvals
    out[[2]] <- phat
    return(out)
}

## Figure 2, left panel
sonfig <- pfig(son4, d)
plotphat <- apply(sonfig[[2]], 2, function(x) quantile(x, c(0.025, 0.5, 0.975)))
uvals <- sonfig[[1]]

pdf("figure2left.pdf", height = 5, width = 5)
par(mar = c(5, 5, 1.2, 1))
plot(uvals, plotphat[2, ], type  = "n", frame = FALSE, xlab = "", ylab = "", ylim = c(0, 0.8), xaxt = "n", yaxt = "n", main = "Son as Successor",
     cex.axis = 1.25, cex.main = 1.5)
lines(uvals, plotphat[2, ])
lines(uvals, plotphat[3, ], lty = "dashed")
lines(uvals, plotphat[1, ], lty = "dashed")
axis(2, at = seq(0, 0.8, by = 0.1), labels = seq(0, 0.8, by = 0.1), las = 2, cex.axis = 1.25)
axis(1, at = seq(0, 4.5, 0.5), labels = round(exp(seq(0, 4.5, 0.5)), digits = 0), cex.axis = 1.25)
mtext("Pr(Son Successor = 1)", 2, at = 0.55, outer = TRUE, cex = 1.5, line = -1.3)
mtext("Tenure (Years)", at = 0.55, 1, outer = TRUE, cex = 1.5, line = -2)
dev.off()
graphics.off()

## Figure 2, right panel
brofig <- pfig(b4, d)
plotphat <- apply(brofig[[2]], 2, function(x) quantile(x, c(0.025, 0.5, 0.975)))
uvals <- brofig[[1]]

pdf("figure2right.pdf", height = 5, width = 5)
par(mar = c(5, 5, 1.2, 1))
plot(uvals, plotphat[2, ], type  = "n", frame = FALSE, xlab = "", ylab = "", ylim = c(0, 0.8), xaxt = "n", yaxt = "n", main = "Brother as Successor",
     cex.axis = 1.25, cex.main = 1.5)
lines(uvals, plotphat[2, ])
lines(uvals, plotphat[3, ], lty = "dashed")
lines(uvals, plotphat[1, ], lty = "dashed")
axis(2, at = seq(0, 0.8, by = 0.1), labels = seq(0, 0.8, by = 0.1), las = 2, cex.axis = 1.25)
axis(1, at = seq(0, 4.5, 0.5), labels = round(exp(seq(0, 4.5, 0.5)), digits = 0), cex.axis = 1.25)
mtext("Pr(Brother Successor = 1)", 2, at = 0.55, outer = TRUE, cex = 1.5, line = -1.3)
mtext("Tenure (Years)", at = 0.55, 1, outer = TRUE, cex = 1.5, line = -2)
dev.off()
graphics.off()

## Figure 3
ndfig <- pfig(nd4, d)
plotphat <- apply(ndfig[[2]], 2, function(x) quantile(x, c(0.025, 0.5, 0.975)))
uvals <- ndfig[[1]]

pdf("figure3.pdf", height = 5, width = 5)
par(mar = c(5, 5, 1.2, 1))
plot(uvals, plotphat[2, ], type  = "n", frame = FALSE, xlab = "", ylab = "", ylim = c(0, 0.25), xaxt = "n", yaxt = "n", main = "Successor Deposed", cex.main = 1.35)
lines(uvals, plotphat[2, ])
lines(uvals, plotphat[3, ], lty = "dashed")
lines(uvals, plotphat[1, ], lty = "dashed")
axis(2, at = seq(0, 0.25, by = 0.05), labels = seq(0, 0.25, by = 0.05), las = 2, cex.axis = 1)
axis(1, at = seq(0, 4.5, 0.5), labels = round(exp(seq(0, 4.5, 0.5)), digits = 0), cex.axis = 1)
mtext("Pr(Successor Deposed = 1)", 2, at = 0.55, outer = TRUE, cex = 1.35, line = -1.3)
mtext("Tenure (Years)", at = 0.55, 1, outer = TRUE, cex = 1.35, line = -2)
dev.off()
graphics.off()

## Figure 4
library(foreign)
GG <- read.dta("full_parliaments.dta")
GG$time <- as.numeric(as.character(GG$time))
GG$presence <- as.numeric(as.character(GG$presence))
GG$lntime <- log(GG$time + 1)
npfit <- loess(presence ~ lntime, data = GG)
nd <- seq(min(GG$lntime), max(GG$lntime), by = 0.1)
names(nd) <- "lntime"
lo_pred <- predict(npfit, newdata = as.matrix(nd), type = "response", se = T)
pdf("figure4.pdf", height = 7, width = 7)
plot(nd, lo_pred$fit, type="l", ylim = c(0.25, 0.5), axes = F, xlab = "Year in Office", ylab = "Probability of a Parliamentary Meeting")
lines(nd, lo_pred$fit + 1.96 * lo_pred$se.fit, lty = 2)
lines(nd, lo_pred$fit - 1.96 * lo_pred$se.fit, lty = 2)
axis(2)
axis(1,at = log(c(1, 2, 3, 4, 7, 12, 30, 20, 33, 55, 72)), labels= c(1, 2, 3, 4, 7, 12, 30, 20, 33, 55, 72))
dev.off()
graphics.off()

## Figure 5
load("parliament_binom.RData")
x <- log(1:73)
prob_predict <- predict(binom1, newdata = data.frame(lnduration = x), type = "response", se.fit = T)
pdf("figure5.pdf", height = 7, width = 7)
plot(log(1:73), prob_predict$fit,lty=1, type = "l", ylab = "Probability of a Parliamentary Meeting",xlab="Tenure (Years)", ylim =c(0.15, 0.6), axes = F)
lines(log(1:73), prob_predict$fit - 1.96 * prob_predict$se.fit, lty = 2)
lines(log(1:73), prob_predict$fit + 1.96 * prob_predict$se.fit, lty = 2)
axis(2)
axis(1, at = log(c(1, 2, 3, 4, 7, 12, 30, 20, 33, 55, 72)), labels= c(1, 2, 3, 4, 7, 12, 30, 20, 33, 55, 72))
dev.off()
graphics.off()
